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Abstract. In this paper wc qualitatively study radial solutions of the somilin- 
ear elliptic equation Au + n" = with u{0) = 1 and u'{0) = on the positive 
real line, called the Emden-Fowler or Lane-Emden equation. This equation is 
of great importance in Newtonian astrophysics and the constant n is called 
the polytropic index. 

By introducing a set of new variables, the Emden-Fowler equation can be 
written as an autonomous system of two ordinary differential equations which 
can be analyzed using linear and nonlinear stability analysis. We perform 
the study of stability by using linear stability analysis, the Jacobi stability 
analysis (Kosambi-Cartan-Chern theory) and the Lyapunov function method. 
Depending on the values of n these different methods yield different results. 
We identify a parameter range for n where all three methods imply stability. 



1. Introduction 

The study of equations of the form Am + f{u) = in bounded or unbounded 
domain has received great attention over the last decades. Many natural processes 
are described by equations of this type. Of particular interest for astrophysics 
are the radial solutions of the Emden-Fowler equation where f{u) = u" and the 
basic equation becomes a second order non-linear ordinary differential equation. 
The mathematical properties have been extensively studied in the past, see for 
instance [IllMlIIlIiaililllllllliaiinillllllllTliaililHlinil In 
particular, the qualitative phase space analysis of the dynamical system associated 
to the Emden-Fowler equation proved to be very useful for the understanding of 
the general mathematical and physical properties of its solutions. 

It is the purpose of the present paper to consider the qualitative study of the 
stability properties of the Emden-Fowler equation which can be transformed to 
an autonomous system of differential equations. After using linear perturbation 
theory to analyze the stability of the critical points, we investigate these equations 
in geometric terms by using the general path-space theory of Kosambi-Cartan-Chern 
(KCC-theory), inspired by the geometry of a Finsler space [28j [TTJ [13]. The key 
idea behind this approach is to geometrize the differential equation itself, rather 
than analyzing its solutions directly. This yields to a manifold (Finsler space) 
similar to phase-space, whose curvature properties determine the behavior of all 
solutions in a non-linear setting. The KCC theory has been applied for the study 
of different physical, biochemical or technical systems (see O [U IH |3l 13 |32l |33l 
|371 mi a-nd references therein). Moreover, to complete our study, we also use the 
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Lyapunov function method to analyze the stabihty properties. By combining these 
three methods, we identify a previously unrecognized parameter range for n where 
stability is implied by all methods. 

The paper is organized as follows. The basic equations and the linear perturba- 
tion theory are presented in Section [2] We review the mathematical formalism of 
the KCC theory and apply it in Section [3l The Lyapunov function method is used 
in Section IH We discuss and conclude our results in Section [S] 

2. Basic equations 

2.1. Emden-Fowler equation. The properties of the static Newtonian astrophys- 
ical objects can be completely described by the gravitational structure equations, 
which are the mass continuity equation and the equation of the hydrostatic equi- 
librium. 

Definition (Gravitational structure equations of Newtonian stars). The basic 
equations describing the equilibrium properties of Newtonian stars are given by 

(2.1) ^^4M^).^ 

ar 

(2.2) = -^Pir), 

ar r-^ 

where p{r) > is the density, p{r) > is the thermodynamical pressure and m(r) 
is the mass inside radius r, respectively, satisfying the condition m(r) > 0,Vr > 0. 

The system of differential equations, given by Eqs. (|2.1[) and (|2.2[) . contains three 
unknown functions, but it consists of only two equations. To close it one must pre- 
scribe an equation of state, p — p{p), which relates the thermodynamical pressure 
to the density of the fluid. We assume that the equation of state is sufficiently 
smooth for all p > {). 

Remark 1. From an astrophysical point of view it seems also well motivated to 
prescribe in some situations the density distribution of matter /o(r), rather than 
the equation of state. However, the general treatment of the problem shows that 
this approach may yield solutions with a non-regular center, which are therefore 
un-physical. 

Definition (Polytropic perfect fluid). An isotropic fluid distribution for which the 
pressure and the density are related by a power law of the form 

(2.3) p = i^pi+i/", 

is called a polytropic perfect fluid. K and n are constants, and ri 7^ is called the 
polytropic index. 

Theorem 2.1 (Emden-Fowler equation). For a polytropic perfect fluid sphere the 
gravitational structure equations V2. 1]) and \2. 2\) are equivalent to the following sin- 
gle second order non-linear differential equation in 0{^), called the Emden-Fowler 
equation 

(2.4) e" + ^e' + 9" ^ 0, 

where the variable 9 is related to the density by means of the definition p — Pcd^ , 
and n is called the polytropic index. The initial conditions for the Emden-Fowler 
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equation are 9{0) = 1 and 9' (0) = 0, respectively, where the prime denotes the 
derivative with respect to the dimensionless radial coordinate ^. 

Proof. By eliminating the mass between the structure equations (|2.1I) and (I2.2p we 
obtain a single second order non-linear differential equation 

. , 1 d fr'^dp\ 

^Tr[jlr)--'''''P- 

It is convenient to represent the density in terms of a new dimensionless variable 
^by 

(2.6) P = Pce", 

giving for the pressure p = Kpl^^^^O^^^^ . Let us furthermore introduce a new 
dimensionless radial coordinate ^ such that 



(2.7) . = of. „=^(!l±Jg^, „^-l. 

In these new variables Eq. (|2.5p takes the form of the Emden-Fowler equation of 
index n. □ 

Remark 2. In the limiting case ti — > 0, the Emden-Fowler equation is solved by 
^iO\n=o ~ 1 ~ For n = 1, the Emden-Fowler equation ()2.4p becomes linear 

and is solved by 9iO\n=i ~ The only known solution to the non- linear 

equation is when n = 5, and it is given by 6'(C)|„=5 -|- 

2.2. The Emden-Fowler equation as an autonomous system. In order to 
study the stability of the equilibrium points of the Emden-Fowler, we rewrite it in 
the form of an autonomous system of differential equations. 

Theorem 2.2. The Emden-Fowler equation {2.4^ is equivalent to the following 
autonomous system of differential equations 

IE - ^' 

(2.9) ^^_2Gi(«;,g), 
with the function G^{'w,g) given by 



(2.10) G\w,q) = ^ 



n - 5 2(3 - ra) „ 
71 — 1 (ri — ly 



Proof. Let us start by introducing a set of new variables {w,t) defined as 

(2.11) 0(e) = se2/(i-"Mo, e = 

where is the value of ^ at the star's surface, and i? > is a constant (which 
generally can be set to one without any loss of generality), the Emden-Fowler 
equation is equivalent with the following second order differential equation, 

d^w n — bdw 2(n — 3) .„ i „ 

In the following we will restrict the range of the polytropic index n to the range 
71 > 1. From a mathematical point of view, Eq. (|2.12p is a second order non-linear 
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differential equation of the form (Pw/dP + 2G^{w,dw/dt) — 0, with given 
by (|2.10p . Note that G^(0, 0) = 0. By introducing a new variable q — dw/dt, 
this second order differential equation is equivalent to the first order autonomous 
system of equations given by Eqs. (|2.8p and (|2.9p . □ 



Remark 3. For the initial Emden- Fowler equation the range of the radial dimen- 
sionless variable ^ is ^ e [0,^s]. The range of the new variable t = ln(^s/^) is 
t G (oo, 0], so that i = at the surface of the star, where £, = £,s, and t oo aX the 
center of the star = 0). Note that if the solution describes an object of infinite 
extend {n > 5), then ^ e [0, oo) and t e (oo, — oo). 

At the surface of the star (the vanishing pressure surface) the function w(i) 
takes the value w{t = 0) = since 9{£,s) = 0, while w'{t = 0) is given by ■w'{t ~ 
0) — 6''(^s)^s^^"^^''/i3. At the center of the star we have limt^oo w (i) — and 
limt_j.oo w' (i) = 0, respectively. 

The critical points of this dynamical system all reside on the g = line. They 
are given by the solutions of the equation 



(2.13) 



2(n- 3) 
"(n- 1)2' 



0. 



Therefore for the critical points Xi — 
firstly find the point 

(2.14) Xo = 

whose value is n independent. 

For w Eq. (|2.13p is equivalent to 

2(n- 3) 



{wo,qo) of the system (|2.8p and (|2.9p we 
= (0,0), 



(2.15) 

Hence, for 1 < n < 3 we find 

(2.16) Xn = 
while for n > 3 we have 

(2.17) Xn = 



{n-iy 



= B" 




2 {n - 3) 
(n- 1)2 



i/(«-i) 



,0 



2.3. Linear stability analysis of the equilibrium points. To characterize the 
nature of the critical points (|2.14p - (l2.17p we use, as a first method, the linear sta- 
bility analysis [lOj . The results of this analysis can be summarized in the following 
theorem. 



Theorem 2.3. Let the dynamical system 112. 8\) and 112. 9\) . with critical points \2.1^^ 
and {2.11^ , be given. Then the point X^ is stable if 1 < n < 3 and unstable 
otherwise. The point A„ is stable for 3 < n < 5 and unstable otherwise. 

Proof. The eigenvalues of the Jacobian matrix of the first derivatives of the dynam- 
ical system. 



(2.18) 





2(n-3)/(?i-l)2-7iS'' 



1 

(n-5)/(n- 1) 
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are given by 

(2.19) X±= ~ ^ , ± - Vl - 4?ig"-iw"-i. 

At the point Xq we therefore find 

(2.20, = 

which gives the following for the eigenvalues 

(2.21) A+ < 0, A_ < 0, 1< n < 3, 

(2.22) A+ > 0, A_ < 0, 3 < n. 

These eigenvalues characterize the critical point Xq, and we have 
1 < n < 3 nodal sink (stable) 
3 < n saddle point (unstable). 

At Xn (here we refer to all n > 1) one finds 

(2.23) A± = — — -(n~5±y/l + 22n-7n^). 

2(n — 1) 

We obtain the following for the eigenvalues 

(2.24) A+ > 0, A_ < 0, 1 < n < 3, 

(2.25) A+<0, A_ < 0, 3 < n < (11 + 8\/2)/7, 

(2.26) ReA+<0, ReA_<0, (11 + 8%/2)/7 < ?i< 5, 

(2.27) ReA+>0, ReA_>0, 5 < n. 
These eigenvalues characterize the critical point Xn, and we find 

1 < n < 3 saddle point (unstable) 

3 < n < (11 + 8V2)/7 nodal sink (stable) 
(11 + 8%/2)/7 < n < 5 spiral sink (stable) 

5 < n spiral source (unstable). 

This completes the linear stability analysis. □ 

3. KOSAMBI-CARTAN-ChERN theory AND THE JACOBI STABILITY 

3.1. Theory. We recall the basics of Kosambi-Cartan-Chern (KCC) theory to be 
used in the sequel [llIlinilllllSllSnilMlISS]- Let 7W be a real, smooth n- 
dimensional manifold and let TAi be its tangent bundle. Let us choose (x*) = 
(x^, cc^, . . . , x"), (y*) = (j/^, 2/^, . . . , J/") and the time i be a 2n + 1 coordinates 
system of an open connected subset $7 of the Euclidian (2n + 1) dimensional space 
R" X R" X R\ where 

' dx^ dx^ c?x"~ 
dt ^ dt ' ^ dt 



(3.1) f 



We assume that t is an absolute invariant, and therefore the only admissible change 
of coordinates will be 

(3.2) t = t, =^'(x^x^...,x"), i G {l,2,...,n}. 
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The equations of motion of a dynamical system can be derived from a Lagrangian 
L via the Euler-Lagrange equations. For a regular Lagrangian the Euler-Lagrange 
equations are equivalent to a system of second-order differential equations 

(3.3) ^ + 2G\x\y\t)^0, z e {1, 2, . . . , n}, 

where each function G* {x^ , , t) is C°° in a neighborhood of some initial conditions 
(xq, yo, to) in D,. The system of Eqs. p.3p is equivalent to a vector field (semi-spray) 
S, where 

(3.4) S = y'— -2GHx\y\t) — , 
which determines a non-linear connection Nj, defined as 



(3.5) N'j ^ 



dyr 

More generally, one can start from an arbitrary system of second-order differential 
equations of the form p.Sp . with no a priori given Lagrangian, and study the 
behavior of its trajectories, by analogy with the trajectories of the Euler-Lagrange 
equations. 

If the coordinate transformations given by Eq. (j3.2p is nonsingular, then the 
KCC-covariant differential of a vector field ^'(x) on the open subset il C M" x M" x 
R-^ is defined as 

(«' + 

For = we obtain 

(3.7) -^=N]y^ -2G^ = ^e\ z € {1, 2, . . . , n} 

The contravariant vector field on Q, is called the first KCC invariant. 

Let us now vary the trajectories x^{t) of the system (|3.3I) into nearby ones, 
according to 

(3.8) x\t)^x\t)+7^C{t), 

where |?7| is a small parameter, and ^*(i) are the components of some contravariant 
vector fields, defined along the path x'^{t). Substituting Eqs. (|3.8|) into Eqs. (|3.3|) . 
and taking the limit 77 — 0, we obtain the variational equations 



The covariant form of Eq. p.9p is 

di2 



(3.10) -i^-PK'. 



where 

dG' , ,dm . , dm 

(3-11) P'; = -^inr ~ ^G'G], + y'^ + NIN ' = 



and 



(3.12) G}, 



dy^ ' 
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Eq. p.lOp is the Jacobi equation, G*; is called the Berwald connection, and Pj is 
the second KCC-invariant, or the deviation curvature tensor O [211 [321 123] ■ When 
the system (I3.3|) describes the geodesic equations in either Riemann or Finsler 
geometry, Eq. (|3.10p is the Jacobi field equation. 

Remark 4. One can also define higher order invariants of the system p.3p . The 
third, fourth and fifth invariants are given by 

1 / dP^ dPl\ dPL dG\^ 

The third invariant is interpreted as a torsion tensor, while the fourth and fifth 
invariants are the Riemann- Christoffel curvature tensor, and the Douglas tensor, 
respectively [T]. In a Berwald space these tensors always exist, and they describe 
the geometrical properties of a system of second-order differential equations. 

The Jacobi stability for a dynamical system is defined as follows [H [351 [33] : 

Definition (Jacobi stability). If the system of differential equations p.3p satisfies 
the initial conditions ||a;*(to) — i*(io)|| = 0, ||a;*(io) — a;H*o)|| 0; "^ith respect to the 
norm || • || induced by a positive definite inner product, then the trajectories of p.3p 
are Jacobi stable if and only if the real parts of the eigenvalues of the deviation 
tensor Pj are strictly negative everywhere, and Jacobi unstable, otherwise. 

3.2. Application to the Emden-Fowler equation. 

Theorem 3.1. Let the dynamical system i2. 8\) and Ii2.9\) . with critical points ^2.14^ 
and ^2.17\ l, be given. Then the point Xq is Jacobi unstable. However, if n > 
(11 + 8\/2)/7, the point Xn is Jacobi stable, and unstable otherwise. 

Proof. By denoting := w and := dx'^/dt = dw/dt, Eq. (|51^ can be written 
as 

(3.14) ^ + 2G\x\y') = 0, 



where 

(3.15) G\x\y' 



2(n-3) 1 



y'-?^^^' + B--'ix'r 



n — l" (n — 1)^ 

respectively, which can now be studied by means of the KCC theory. As a first step 
in the KCC stability analysis of the Newtonian polytropes we obtain the nonlinear 
connections = dG^ /dy^, associated to Eqs. (|3.14l) . and which is given by 

(3.16) Nl{x\y') = -l^. 

2 n — 1 

For the Emden-Fowler equation the associated Berwald connection identically van- 
ishes 

(3.17) = ^ = 

Finally, the second KCC invariant, or the deviation curvature tensor P^, defined 

as 
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is given by 

(3.19) Pl{x\y') = \-nB''^-\x'r~\ 
In the initial variables is given by 

(3.20) Pl{w,q) = i - 

Evaluating Pl{w,q) at the critical points X„, given by Eqs. (|2.14l) - (|2.17l) . we 
obtain 

(3.21) pi(Xo) = i>0, Vn, 

Therefore, according to the Jacobi stability theorem, the point Xq is unstable 
since Pi > 0. However, for n > (11 + 8^/2)/ 7 the function becomes strictly 
negative and Jacobi stability is established. □ 

3.3. Physical interpretation of Jacobi stability. In the initial variables the 
deviation curvature tensor P^ is given by 

(3.23) pi(^,e)= 1 

P^ can be expressed in a simple form with the use of Milne's homological variables 
{u,v), see [mils], defined as 

with the use of whom we obtain 

(3.25) Pl{u,v) ^^~nuv. 

From a physical point of view u(r), defined as u(r) = dlnTO(r)/dlnr = 3p(r)/ p(r), 
is equal to three times the density of the star at point r, divided by the mean den- 
sity of matter, p{r), contained in the sphere of radius r. The variable w(r), defined 
as v{r) = -dliip{r)/dliir = {3/2)[Gm{r)/r]/[3p{r)/2p{r)], is 3/2 of the ratio of 
the absolute value of the gravitational potential energy, \Eg\ — Gm{r)/r, and the 
internal energy per unit mass, Ei = (3/2)p/p, so that v{r) — {3/2)\Eg\/ Ei [23l [8]. 
In terms of these physical variables the deviation curvature tensor is given by 

^^■^^> 2 pir) E^ 

The condition of the Jacobi stability of the trajectories of the dynamical system 
describing a Newtonian polytropic star, P^ < 0, can therefore be formulated as 



(3.27) 
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4. LYAPUNOV FUNCTION STABILITY ANALYSIS 



4.1. Theory. When discussing the mathematical background, we closely follow [35l 



Deflnition (Lyapunov function). Given a smooth dynamical system x — f{x), 
X e M", and an equilibrium point xq, a continuous function V : M" — ;> i? in a 
neighborhood U of xq is a Lyapunov function for the point xq if 

(1) V is differentiable in t/ \ {a;o} 

(2) V{x) > Vixo) 

(3) V{x) < for every x(EU\ {xq}. 

The neighborhood U is constrained by the third condition. When 3. holds, 
the method provides information not only about the asymptotic stability of the 
equilibrium point, but also about its basin of attraction, which must contain the 
set U. 

The existence of a Lyapunov function V for an autonomous system of differential 
equations guarantees the stability of the point xq. This qualitative result can be 
obtained without explicitly solving the equations, and this information cannot be 
obtained by using linear stability analysis or first order perturbation theory. The 
following theorem holds: 

Theorem 4.1 (Lyapunov stability). Let Xq be an equilibrium point of the system 
X = fix), where f : U ^ M" is locally Lipschitz, and U C M" is a domain that 
contains xq. If V is a Lyapunov function, then 

(1) if V{x) = dV/dxf is negative semi- definite, then x — xq is a stable equi- 
librium point, 

(2) if V{x) — dV/dxf is negative definite, then x ~ Xq is an asymptotically 
stable equilibrium point. 

Moreover, if \\x\\ — >■ oo implies that V{x) — >■ oo for all x, then xq is globally stable 
or globally asymptotically stable, respectively. 

If the condition 3. of the Lyapunov function definition holds strictly, the existence 
of a Lyapunov function has important consequences for the behavior of the time 
dependent perturbations of autonomous systems. 

4.2. Application to Emden-Fowler equation. 

Theorem 4.2. Let the dynamical system k2. and i2. 9\) . with critical points {2.1^^ 
and {2.11^ , be given. The equilibrium state is asymptotically stable is 1 < n < 3, 
and the equilibrium state Xn is asymptotically stable is 3 < n < b. 

Proof. One possible Lyapunov function V{w, q) associated to the system given by 
Eqs. I|2.8p and (|2.9p can be chosen following [S^. Following the variable gradient 
method (this means setting WV — f G M"), we set 



such that the critical points correspond to WV — (since this corresponds to / = 0). 
This yields the following Lyapunov function 



(4.1) 




{n-iy 



B 



m — 1 



(4.2) 



V{w,q) = -q^ 




n + 1 



w 



10 



C. G. BOHMER AND T. HARKO 



By definition, the Lyapunov function must have a local niinimum at the critical 
points. To check this, we consider the Hessian of (|4.2|) . which is given by 

(4.3) ^(^)^(^-2(n-3)/(n-l)^ + 5"-n."- 0^ 
and has eigenvalues 

(4.4) Ai = -2(n-3)/(n- l)2 + B"-inw"-\ A2 = 1 
At the point Xq — (0, 0) for n > 1 we have 

(4.5) Ai = -2(n-3)/(n-l)2, A2 = 1, 

and hence there is a local minimum near Xq if 1 < n < 3. 
At the point Xn, on the other hand, we find 

(4.6) Ai =2(n-3)/(n-l), A2 = 1, 

and hence there is a local minimum near Xn if n > 3, and this function has no local 
minimum if the index satisfies 1 < rt < 3. 
The Lyapunov function (|4.2|) satisfies 

dV dV dw ^ dV dq n — 5 2 
dt dw dt dq dt n — 1 ' 
and therefore we find that 

(4.8) V <0, 1 < n < 5. 

Hence, according to the Lyapunov theorem, the equilibrium state Xq is asymptot- 
ically stable equilibrium points for 1 < n < 3 while the equilibrium state Xn is 
asymptotically stable equilibrium points for 3 < n < 5. Global asymptotic stability 
results cannot be inferred from our Lyapunov function. □ 

Remark 5. Now we are faced with the difficulty of finding a Lyapunov function for 
other values of n. Recall that the Lyapunov theorem is in fact a quite restrictive 
statement, though very powerful. If stability properties cannot be established us- 
ing a particular Lyapunov candidate function, it does not mean that the point is 
unstable. 

5. Discussions and final remarks 

We considered the stability properties of the Emden-Fowler equation. For the 
stability analysis we used three methods: linear stability analysis, the so-called 
Jacobi stability analysis, or the KCC theory, and the Lyapunov function method. 
Our results show that the study of the stability properties of the Emden-Fowler 
equation turns out to be a mathematically interesting problem. 

The study of the stability has been done by analyzing the behavior of the steady 
states Xq and X„. It is most interesting to note that these three method for some 
parameter region of the index n give different stability properties. Our results can 
be summarized in Table [1] 

Hence, we have found that a solution of the Emden-Fowler equation can be 
regarded as stable if and only if its polytropic index satisfies (11 -l-Sv^)/? < n < 5. 
For other parameter ranges the various methods used yield different results. It is 
however possible to give a good geometrical interpretation of the situation 3 < n < 
(11 -|- 8V2)/7 when the linear perturbation theory and the Lyapunov method both 
imply stability while the point is Jacobi unstable. 
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index n/method 


Linear 


Jacobi 


Lyapunov function 


1 < n < 3 


unstable 


unstable 


inconclusive 


3 < ?i< (11 + 8V2)/7 


stable 


unstable 


stable 


(ll + 8\/2)/7<n<5 


stable 


stable 


stable 


5 < n 


unstable 


stable 


inconclusive 



Table 1. Stability properties of polytropic perfect fluid spheres 
with polytropic index n as implied by linear perturbation theory, 
the Jacobi stability method (KCC theory), and the Lyapunov func- 
tion method. 



Let us recall [32] that the Jacobi stability of a dynamical system is regarded 
as the robustness of the system to small perturbations of the whole trajectory. 
This is a very convenient way of regarding the resistance of limit cycles to small 
perturbation of trajectories. It is interesting to remark that the stable nodal sink 
(3 < n < (ll-|-8-\/2) /7) in linear perturbation theory) is actually Jacobi unstable. In 
other words, even though the system's trajectories are attracted by the critical point 
Xn one has to be aware of the fact that they are not stable to small perturbation 
of the whole trajectory. 

We conclude by mentioning that our method can in principle be applied to all 
radially symmetric problems of the form Au + f{u) = 0. 

Acknowledgments. We thank Steven Bishop for the useful discussions and valu- 
able comments on the manuscript. The work of T. H. was supported by an RGC 
grant of the government of the Hong Kong SAR. 
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